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ABSTRACT 

We investigate the power of wavelet techniques in detecting non-Gaussianity in the 
cosmic microwave background (CMB). We use the method to discriminate between 
an inflationary and a cosmic strings model using small simulated patches of the sky. 
We show the importance of the choice of a good test statistic in order to optimise the 
discriminating power of the wavelet technique. In particular, we construct the Fisher 
discriminant function, which combines all the information available in the different 
wavelet scales. We also compare the performance of different decomposition schemes 
and wavelet bases. For our case, we find that the Mallat and a trous algorithms are 
superior to the 2D-tensor wavelets. Using this technique, the inflationary and strings 
models are clearly distinguished even in the presence of a superposed Gaussian com- 
ponent with twice the rms amplitude of the original cosmic string map. 

Key words: methods: data analysis-cosmic microwave background. 
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INTRODUCTION 

The Cosmic Microwave Background (CMB) is currently one 
of the most powerful tools of cosmology to investigate the 
theories of structure formation in the Universe. In particular, 
the standard inflationary model predicts Gaussian fluctua- 
tions whereas alternative theories such as topological defects 
scenarios give rise to non-Gaussianity in the temperature 
field. Therefore, the study of the Gaussian character of the 
CMB will allow us to discriminate between these two com- 
petitive theories of structure formation. In addition, future 
high-resolution maps should be investigated for any trace of 
non-Gaussianity that may be imprinted by foregrounds or 
systematics. 

In order to test the Gaussianity of the CMB, a large 
number of methods have been proposed (for a review see 
Barreiro 2000). Perhaps the most straightforward of these 
is the measurement of the skewness and kurtosis of the dis- 
tribution of the pixel temperatures in the map (Scaramella 
& Vittorio 1991). More complex methods include proper- 
ties of hot and cold spots (Coles & Barrow 1987; Martinez- 
Gonzalez et al. 2000; Barreiro, Martinez-Gonzalez & Sanz 
2001), extrema correlation function (Naselsky & Novikov 
1995, Kogut et al. 1996, Barreiro et al. 1998, Heavens & 
Sheth 1999), Minkowski functional (Coles 1988, Gott et al. 
1990, Kogut et al. 1996), bispectrum analysis (Ferreira et 
al. 1998, Heavens 1998, Magueijo 2000), multifractals (Pom- 
pilio et al. 1995) and partition function (Diego et al. 1999, 
Martinez-Gonzalez et al. 2000). 

In addition, tests of Gaussianity based on wavelets have 



become quite popular in the last years. Ferreira, Magueijo 
& Silk (1997) propose a set of statistics using cumulants 
and defined in wavelet space. Hobson, Jones & Lasenby 
(1998) investigate the power of the cumulants of the dis- 
tribution of wavelet coefficients at each scale on detecting 
the Kaiser-Stebbins effect using CMB simulations of small 
patches of the sky. A similar method was introduced by Forni 
& Aghanim (1999) and tested using simulated CMB maps 
containing secondary anisotropics (Aghanim & Forni 1999) 
and on the COBE data (Aghanim, Forni & Bouchet 2001). 
Pando, Valls-Gabaud & Fang (1998) calculated the skew- 
ness, kurtosis and scale-scale correlation coefficients of the 
COBE data in the QuadCube pixelisation using planar or- 
thogonal wavelets. This study was extended by Mukherjee 
et al (2000). Barreiro et al (2000) applied the same method 
to the COBE data in HEALPix pixelisation using spherical 
wavelets. 

This paper is an extension of the work done by Hobson 
et al (1998) (HJL, hereinafter). As in HJL, the fourth cumu- 
lant of the distribution of wavelet coefficients is calculated at 
each wavelet scale using high-resolution CMB simulations of 
small patches of the sky for Gaussian and non-Gaussian sim- 
ulations. However, in this work, all this information is com- 
bined into the Fisher discriminant function in order to get a 
more statistically significant result. In addition, we compare 
the performance of three different ways of constructing 2- 
dimensional wavelets: tensor (as in HJL), Mallat and the a 
trous algorithm. Moreover, some of the studied wavelets are 
invariant under rotations of 90, 180 and 270 degrees of the 
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original data, which solves the earlier problem of orientation 
sensitivity. 

The outline of the paper is as follows. gives a 
brief introduction to the three different constructions of 2- 
dimensional wavelets. In §^ our non-Gaussianity test is ex- 
plained in detail, also pointing out the interest of the Fisher 
discriminant function. The results of applying this test to 
simulated data are given in §^ Finally, our conclusions are 
presented in 



2 THE DISCRETE WAVELET TRANSFORM 

The wavelet transform has been extensively discussed else- 
where (e.g. Daubechies 1992). In particular, a clear intro- 
duction to wavelets is given by Burrus, Gopinath & Guo 
(1998). We therefore give only a brief introduction, pointing 
out the differences between the three considered decompo- 
sition schemes. The basis of the wavelet transform is to de- 
compose the considered signal (the temperature field in our 
case) into a set of wavelet coefficients. Each of them con- 
tain information about the signal at a given position and 
scale. Therefore, wavelets allow one to study the structure 
of an image at different scales and at the same time keep lo- 
calisation. There is not a unique choice for a wavelet basis. 
Moreover, different algorithms schemes can be used to con- 
struct 2-dimensional wavelets. In particular, we consider the 
tensor, Mallat (or multiresolution) and a trous algorithms. 

2.1 Tensor algorithm 

This is the type of wavelets used in the analysis of HJL, 
where a more detailed description can be found. In 1- 
dimension, the wavelet basis is constructed from dilations 
and translations of the mother (or analysing) wavelet func- 
tion ip and a second related function called the father (or 
scaling) wavelet function: 

= 2^tP(2 3 - j x~1) , 

4> jlt = 2^ (^"'as - l) , (1) 

where j ^ J — 1 and ^ I ^ 2 J — 1 are integer denoting 
the dilation and translation indices, respectively, and 2 ,J is 
the number of pixels of the considered discrete signal f(x). 

In order to construct a real, orthogonal and compactly- 
supported wavelet basis, such as the Haar and Daubechies 
4 wavelets used in this work, these functions must satisfy 
some highly-restrictive mathematical relations (Daubechies 
1992). The discrete signal f(x) can thus be written as 

f(xi) = oo,o</>o,o(a;i) + ^ X] d i.i^jA x i) > ( 2 ) 
j I 

being a, d the approximation and detail wavelet coefficients 
respectively. These coefficients can be found in a recursive 
way starting from the data vector f(xi) = aj t i , I — 0, ..,2 J — 1 

aj-ij = y^h(m — 2l)aj >m 

m 

dj-x,i = y^g(m - 2l)a j>m . (3) 

m 

where h, g are the low and high-pass filters associated to 



the scaling and wavelet functions through the refinement 
equation (see e.g. Burrus et al. 1998). At each iteration, 
the vector of length 2 J is split into 2 J ~ detail components 
and 2-'~ 1 smoothed components. The decimated smoothed 
components are then used as input for the next iteration 
to construct the detail coefficients at the next larger scale. 
This produces a number of wavelet coefficients equal to the 
original number of pixels. As the index j increases from to 
J — 1, the wavelet coefficients correspond to the structure of 
the function on increasingly smaller scales, with each scale 
a factor of 2 finer than the previous one. 

The extension of the DWT to 2-dimensional signals is 
obtained by taking tensor products of the one-dimension 
wavelet basis as follows 

<fo,0;h,h{x,y) = <t>o,h(x)(t>o,i 2 (y) 
(ji,0\h.h(x,y) = i>ji,li(%)<t>o : l 2 {y) 

£o,j 2 ;h,h{ x ,y) = <t>o,i 1 { x )^ 02 ,i 2 {y) 
i>iiMh,iA x ,y) = Tpii,h( x )i>j2,i2(y)- 

If the two-dimensional pixelised image has dimensions 2 J x 
2 J then we have < ji < J - 1 and < h < 2 J1 - 1, and 
similarly for 32 and fa. 

With this wavelet decomposition scheme, the wavelet 
transform is performed independently in both directions. 
By analogy with the one-dimensional case, as j\ (32) in- 
creases the wavelet coefficients contain information about 
the horizontal (vertical) structure in the image on increas- 
ingly smaller scales. In particular, the scale corresponding to 
each region is approximately ~pixel sizex2 J_Ji in each di- 
rection. Therefore, regions with j\ = 32 contain coefficients 
of two-dimensional wavelets that represent the image at the 
same scale in the horizontal and vertical directions. Con- 
versely, domains with j\ 7^ 32 contain wavelet coefficients 
describing the image on different scales in the two direc- 
tions. 

2.2 Mallat algorithm 

For the previous wavelets, the wavelet transform is applied 
separately for each direction of the image. This means that 
some of the domains contain mixing of structure with differ- 
ent scales ji and 32 in the horizontal and vertical directions. 
However, it is possible to use one-dimensional bases to con- 
struct two-dimensional bases that do not mix scales and can 
be described in terms of a single scale index j (Mallat 1989). 
At scale 3, these bases are given by 

$r,h,i2( x ,y) = (p^iA^^^hiy) 

Tpr,h,i 2 ( x ,y) = 0i,h( x )i'j,h(y) 
ip?;h,h{ x >y) = 4 , j,h( x ) i Pj,h(y)- 

In this wavelet decomposition, the temperature map (scale 
J) is decomposed into an approximation image (that is ba- 
sically a lower resolution version of the original image) plus 
three details images. These detail coefficients contain hori- 
zontal, vertical and diagonal structure at scale J — 1 that 
encode the differences between the original and lower reso- 
lution images. The same decomposition is applied then to 
the approximation image, generating a lower resolution im- 
age and another three details at scale J — 2. This process 
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is repeated down to the lowest resolution level considered. 
The scale of the structure contained by a wavelet coefficient 
with index j correspond approximately to pixel sizex2 / ~-' 
As for the tensor case, the wavelet basis is non-redundant 
and orthogonal. Thus, the number of wavelet coefficients is 
the same as the number of pixels. We also note that the 
diagonal wavelet coefficients obtained with the Mallat al- 
gorithm are the same as those with ji — j% in the tensor 
decomposition. 

Although different orthogonal wavelet bases were tested 
for the tensor and Mallat algorithms, we present our re- 
sults for Haar and Daubechies 4 (see §Q). We note that it 
is not possible to construct a real orthogonal wavelet ba- 
sis with compact support whose wavelet basis functions are 
also symmetric. Therefore the wavelet transform will not be 
invariant under rotations. However, the Haar wavelet func- 
tions are antisymmetric. This means that our test will be 
invariant under rotations of 90, 180 and 270 degrees of the 
data with respect to the original orientation when using the 
Haar wavelets basis. 

2.3 A trous algorithm 

The a trous ('with holes') algorithm (Holschneider & 
Tchamitchian 1990; Shensa 1992) is a redundant non- 
orthogonal wavelet transform. Given a data vector f(xi) = 
cj,(, 1 = 0, 2 J — 1, the approximation coefficients for the 
next lower resolution are obtained as 



i + 2( J -i)m 



(4) 



i.e., each step is a convolution of the data with the filter 
h and varying step sizes 2 J ~K The function ft is a discrete 
low pass filter associated to a scaling function 4>{x) that 
can be derived from a refinement relation (see e.g. Starck 
& Murtagh 1994). The corresponding detail coefficients are 
simply obtained by subtracting the smoother image from 
the original one: 

dj—i,i = Oj.fc — Oj_i,fc, (5) 

Since no decimation is carried out between consecutive filter 
steps, the number of wavelet coefficients at each scale is 2 J . 
Note that the reconstruction of the signal is particularly 
simple in this case, obtained by adding the detail wavelet 
coefficients at all scales plus the approximation coefficients 
at the lowest considered resolution: 
j-i 

CJJ = Co,! + WjJ (6) 
3=0 

We have chosen a £>3-spline for the scaling function, 
which leads to a convolution with a mask (^,i,|,i,-L)in 
1-dimension and 
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(7) 



in 2-dimensions (e.g. Starck & Murtagh 1994). Note that, 
due to the symmetry of the mask, the a trous wavelet trans- 
form is invariant under rotations of the data of 90, 180 and 
270 degrees with respect to the original orientation. 



3 THE NON-GAUSSIANITY TEST 

The non-Gaussianity test used in this paper is basically the 
same as that of HJL. The main idea is to perform a wavelet 
transform of the (possibly) non-Gaussian image and obtain 
certain statistic of the wavelet coefficients in each region 
of the transform. In particular, we estimate the fourth cu- 
mulant K4 for each scale. This quantity is calculated using 
fc— statistics as explained in HJL. We repeat this test for 
a large number (10000) of non-Gaussian and Gaussian re- 
alisations and then compare the distribution of the fourth 
cumulant spectra «4(r) obtained from the two different pop- 
ulations. This is different to HJL where only a single non- 
Gaussian map was used to test the power of the wavelet 
technique (see §^j). Each Gaussian map is obtained from one 
of the non-Gaussian maps by randomising its phases. There- 
fore, Gaussian and non-Gaussian maps share the same power 
spectrum by construction. This is important to ensure that 
any possible differences found between both populations are 
due to the non-Gaussian character of the test map rather 
than differences in the power spectrum. 

As discussed in HJL, due to edge effects, not all wavelet 
coefficients can be used in the analysis. We need to discard 
those coefficients obtained from wavelet functions that cross 
the boundary of the image to avoid misidentifying these 
borders as non-Gaussian features. To identify the affected 
wavelet coefficients, we perform a wavelet transform of a 
map consisting of zeroes except for a border of non-zero 
pixels on the map boundary. The resulting non-null wavelet 
coefficients correspond to edge-crossing basis functions and 
are therefore ignored in the subsequent analysis. 

One can, of course, repeat the test using cumulants 
other than K4. We have obtained results for K3, but find it 
to be a less discriminating statistic than in the case studied 
here. 



3.1 The discriminating power of a test 

In HJL, the significance of the test was given by the level 
of the highest detection of non-Gaussianity. However, a sin- 
gle detection of non-Gaussianity may not be statistically 
significant if all the other estimators are consistent with a 
Gaussian distribution. Conversely, when the level of con- 
fusion is high, it may happen that no single statistic de- 
tects non-Gaussianity at a high significance level. Neverthe- 
less, if we use all the available statistics together, we may 
find a systematic displacement with respect to the Gaussian 
distribution, which would allow us to discriminate between 
both models. Therefore, a more statistically rigorous analy- 
sis should take into account all the information provided by 
the different regions of the wavelet transform. In order to 
do this, we need to construct a test statistic that maximises 
the discriminating power between the two models. 

Suppose that given a set of n measurements x = 
(xi, x n ), we want to check the validity of a hypothesis 
Ho, usually called the null hypothesis, against an alternative 
hypothesis Hi . Each hypothesis specifies a joint probability 
/(x|iJo) and /(x|_Hi) for obtaining the particular values of 
our data. In our case, the measurements correspond to the 
value of the fourth cumulant at each wavelet scale and we 
want to test if our sample is drawn from a Gaussian model or 
from a cosmic strings scenario. In order to compare the level 
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Figure 1. An example of the discriminating power of a test be- 
tween two hypotheses Ho and Hi. The solid line corresponds to 
the probability g(t\Ho) of obtaining a certain value of t if Hq is 
true, whereas the dashed line gives the same probability g(t\Hi) 
if Hi is the correct hypothesis. 



of agreement between the measured data and the considered 
hypotheses, it is useful to construct a function of the mea- 
sured variables called a test statistic £(x). In principle, this 
statistic can be a multidimensional vector, but we consider 
only the case when t is a scalar, since this greatly simpli- 
fies the problem by reducing the number of dimensions of 
the data whilst keeping enough information to discriminate 
between the considered hypotheses. The probability of ob- 
taining a certain value of t will be given by g(t\Ho) if Ho is 
true and by g(t\Hi) if the correct hypothesis is Hi. In order 
to reject or accept the hypothesis Ho, one can define a crit- 
ical value for t (see Fig. hi). The critical value t c is defined 
so that the probability of observing t under the assumption 
of the hypothesis Ho is less than a certain value a, which 
is called the significance level of the test. In the example of 
Fig. |l] this would be given by 



g{t\Ho)dt 



(8) 



Thus, there is a probability a of rejecting Hq when it is 
actually true. This is called an error of the first kind. An 
error of the second kind occurs when t falls in the accep- 
tance region but it was actually drawn from an alternative 
hypothesis Hi. In our example, in Fig. [l], the probability of 
such an error happening is given by 



g(t\Hi)dt 



(9) 



The quantity p = 1 — /3 is called the power of the test to 
discriminate against the alternative hypothesis Hi at the 
considered significance level. In our wavelet analysis, we will 
compare the values of p for a significance level of a = 0.01 
for each of the studied cases. 



3.2 The 



X 



test 



The key issue now is the choice of a test statistic. One possi- 
bility is to uncorrelate the variables and then obtain a stan- 
dard \ 2 statistic. The procedure is as follows. First, the em- 
pirical covariance matrix V of the variables is obtained from 



a large number of Gaussian simulations. Then, a new set of 
variables x is constructed by a rotation: 



5 T x 



(10) 



where S is a matrix whose columns are the normalised eigen- 
vectors of the covariance matrix V. The same rotation is 
then applied to the non-Gaussian variables. Finally, the \ 2 
statistic is obtained for each Gaussian and non-Gaussian re- 
alisation by calculating 



E 



(ii) 



where x\ corresponds to the ith (Gaussian or non-Gaussian) 
rotated variable, and a[ and \i i to the dispersion and mean 
value of the ith uncorrelated variable for the Gaussian dis- 
tribution. 

This method is actually equivalent to perform a gen- 
eralised x 2 (such as that in Barreiro, Martinez-Gonzalez & 
Sanz 2001) from the original correlated variables that takes 
into account the correlations between them. However, our 
procedure has some advantages. By uncorrelating the Gaus- 
sian variables, we can have some insight into the significance 
of the single detections of non-Gaussianity. In addition, it 
gives us the possibility of investigating a further test statistic 
as follows. We can construct the probability density function 
P(x'i) for each of the uncorrelated variables obtained from a 
large number of Gaussian simulations. Although the require- 
ment that the variables are uncorrelated does not imply that 
they are independent (i.e. that their joint probability is sep- 
arable), we may still construct the test statistic 



(12) 



If the joint probability of x' is a multidimensional Gaus- 
sian, rj 1 is equivalent to the \ 2 (they coincide except for 
an additive constant and a scaling factor). The construction 
of a non-Gaussian \ 2 was already introduced by Ferreira, 
Magueijo & Gorksi (1998) to study the normalised bispec- 
trum of the COBE data. Their definition includes some extra 
factors to ensure that rj 2 defaults exactly to the standard \ 2 
when the considered variables follow a multivariate Gaus- 
sian distribution. 

Unfortunately, we could not estimate this quantity for 
our non-Gaussian realisations. This is due to the fact that 
some of the values of the variables in the non-Gaussian case 
lie far beyond a few cr's with respect to the distribution of the 
Gaussian population. This makes it very difficult to evaluate 
the quantity rf for the non-Gaussian distribution. Instead, 
we compute the standard \ 2 statistic (|l"l| ) and compare its 
performance to the test described below. 



3.3 The Fisher discriminant function 

In the previous tests, we are comparing the level of agree- 
ment of a set of data with some theoretical model. Therefore 
only information about this model is included in the test. 
However, if we believe that our data are most likely drawn 
either from a hypothesis Ho or from an alternative one Hi, 
it is well known that an optimal test statistic in the sense of 
maximum power for a given significance level is given by 
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t(x) = 



/(x|gp) 
/(x|ffi) 



(13) 



which is called the likelihood ratio (Brandt 1992). Unfortu- 
nately, in order to construct this statistic, we need to know 
the joint probabilities /(x|_ffo) and /(x|_H"i), which are in 
general not available. 

A simpler possibility is to construct a test statistic using 
the Fisher linear discriminant function (Fisher 1936). This 
is the optimal linear function of the measured variables for 
separating the probabilities g(t\Ho) and g(t\Hi) in the sense 
of maximum distance between the expected mean values of 
t for each model and minimum dispersion of each of them. 
The test statistic is given by 



t(x) 



with 



a cx W 1 (/j,o - Mi) 
r U = (Vo + Vi)ij 



(14) 



(15) 



where fik and Vk are the expected mean value and covariance 
matrix of the measured data if drawn from the hypothesis 
Hk. Thus i(x) is determined up to an additive constant and 
a scaling factor. 

Actually if /(x|i?o) and f(x\Hi) are both multivari- 
ate Gaussian with common covariance matrices, the Fisher 
discriminant function coincides with the likelihood ratio 
(Cowan 1998). 

In our case, the set of variables x are the fourth cu- 
mulant of the wavelet coefficients at each scale, whereas the 
hypotheses Ho and Hi correspond to our simulated data be- 
ing drawn from a Gaussian and a non-Gaussian population 
respectively (see sQ). In order to estimate the power of our 
non-Gaussianity test, we construct and compare the Fisher 
discriminant function t for each Gaussian and non-Gaussian 
realisation. The mean values and covariance matrices of our 
variables for both models needed to generate the W matrix 
are estimated from a large number of simulations. 



4 APPLICATION TO SIMULATED DATA 

In this section, we apply the wavelet analysis described 
above to the detection of non-Gaussianity in simulated CMB 
maps. As our test map, we use the same realisation of CMB 
anisotropies as in HJL. This is a simulated CMB map due 
to the Kaiser-Stebbins effect from cosmic strings (Kaiser & 
Stebbins 1984). The size of the simulation is 6.4 x 6.4 deg 2 
with a total of 256x256 pixels of size 1.5 arcminutes. 

In order to test the power of our non-Gaussianity test, 
we add to the cosmic strings map a Gaussian component 
that shares the same power spectrum as the test map. The 
weight of the Gaussian component is controlled by a param- 
eter b, that gives the ratio of the non-Gaussian to Gaussian 
proportions as 1:6. We investigate which is the maximum 
level of this Gaussian component (i.e., the largest value of 
b) that can be introduced and still detect the underlying 
non- Gaussianity. 

To make more realistic simulations, the test map is 
convolved with a Gaussian beam of FHWM=5 arcminutes. 
Finally, Gaussian white noise is added with an rms level 



equal to one-tenth that of the convolved map. This de- 
gree of smoothing and level of noise is typical of what may 
be achieved from CMB observations by the Planck satel- 
lite of ESA after foreground removal is performed using the 
Maximum-Entropy method (Hobson et al 1998). 

As pointed out in the previous section, a single non- 
Gaussian map was used to test the power of the wavelet 
analysis in HJL. However, even for the same underlying non- 
Gaussian map, the particular value of our statistics depends 
on the specific Gaussian component and noise realisation 
that are being used. Therefore, the conclusions derived from 
our test will vary from one Gaussian realisation to other. To 
take this into account, we have analysed a large number 
(10000) of non-Gaussian test maps that contain the same 
cosmic strings map but different Gaussian component and 
noise realisations. Ideally, each non-Gaussian map would 
also contain a different cosmic strings realisation. Unfortu- 
nately, strings maps are difficult to simulate and very com- 
putationally expensive. Therefore, a large number of stringy 
maps are not currently available to us and we account only 
for the dispersion introduced by the Gaussian components 
of the map. Nonetheless, although the particular level of 
discrimination achieved by the test may vary for a differ- 
ent choice of a cosmic strings map, we expect the general 
conclusions of this work to remain valid. 

Another important point is the choice of the type of 
wavelet and the wavelet basis. We compare the performance 
of three different ways of constructing wavelets: tensor, Mal- 
lat and a trous algorithm. For the tensor and Mallat algo- 
rithms, we also consider different wavelet bases. In HJL, 
an entropy criterion was used for finding the wavelet basis 
that gives the single highest detection of non-Gaussianity 
at a particular scale. In this paper, however, we perform a 
joint analysis of the information provided by all the scales 
and therefore we cannot apply the same criterion. We have 
tried Haar, Daubechies 4,6,12 and 20, Coifiet 2 and 3 and 
Symmlet 6 and 8 wavelet bases and find that, using the full 
analysis, most of them lead to similar results. In particular, 
more complex wavelets (such as Daubechies 12 and 20) tend 
to discriminate between the models slightly worse, probably 
because of a larger number of wavelet coefficients being dis- 
carded due to edge effects. Therefore, we present our results 
for Daubechies 4 and Haar wavelets since these are two of 
the simplest and also best known orthogonal wavelet bases. 
Moreover, the Haar basis functions are antisymmetric and 
the results of our statistic are therefore invariant under ro- 
tations of the data (by 90, 180 and 270 degrees). 

4.1 Tensor algorithm 

2D-tensor wavelets are constructed in the manner explained 



in §2.1. These are the type of wavelets used in the origi- 
nal work of HJL. In this section we present the results for 
this type of wavelets for the Daubechies 4 and Haar wavelet 
bases. We also show the discriminating power of the Fisher 
discriminant function versus a \ 2 analysis and address the 
orientation sensitivity problem. 

4-1.1 Daubechies 4 

The top panel of Fig. |^ shows the fourth cumulant spectrum 
for the test map with non-Gaussian:Gaussian proportions 
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Figure 2. Top panel shows the fourth cumulant spectra obtained 
from 10000 Gaussian and non-Gaussian simulations for b = 0, us- 
ing tensor Daubechies 4 wavelets. The corresponding uncorrelated 
variables are given in the middle panel. The bottom panel shows 
the histogram of the Fisher discriminant obtained from the Gaus- 
sian (solid line) and non-Gaussian (dashed line) simulations (see 
text for details). 



(1:0) using Daubechies 4 and the original orientation of the 
data. Solid squares and open circles correspond to the non- 
Gaussian and Gaussian values respectively. The Gaussian 
values have been normalised to zero mean value and unit 
dispersion to allow for a straightforward reading in units 
of the Gaussian dispersion of the distance between the two 
models. The error bars show the 68, 95 and 99 per cent 
confidence levels of the Gaussian and non- Gaussian distri- 
butions derived from 10000 simulations in each case. The 
scale corresponding to the region number is given in Ta- 
ble [ij a higher region number corresponds, in general, to 
smaller scales. We see that there are several large detections 
of non-Gaussianity, the highest single detection being ~ 60a 



Table 1. Correspondence of region number with scale for tensor 
wavelets. The scale is calculated as 1.5'(2 J— J1 X 2 J ~ J2 ), where 
2 J X 2 J is the total number of pixels in the map. The fifth and 
sixth columns give the total number of coefficients in the region 
and how many of those coefficients are used to calculate the fourth 
cumulant of their distribution when using Daubechies 4. For the 
Haar wavelet, there is not need of correction for boundary effects 
and all coefficients are used. Finally, the last column indicates the 
number of the corresponding region for Mallat wavelets, for those 
regions that are common for both types of wavelets. For sym- 
metry, wavelet coefficients in domains with j'2 > ji are included 
in the same region as those with scales (j2,ii). Following HJL, 
wavelet coefficients with ji = are not included in the analysis. 



region 


h 


h 


~scale(') 


no.coeff. 


used (D4) 


Mallat 


1 


l 


1 


192 X 192 


4 





3 


2 


2 


1 


96 x 192 


16 







3 


2 


2 


96 x 96 


16 





6 


4 


3 


1 


48 x 192 


32 







5 


3 


2 


48 x 96 


64 


10 




6 


3 


3 


48 x 48 


64 


25 


9 


7 


4 


1 


24 x 192 


64 







8 


4 


2 


24 x 96 


128 


26 




9 


4 


3 


24 X 48 


256 


130 




10 


1 


4 


24 X 24 


256 


169 


12 


11 


5 


1 


12 X 192 


128 







12 


5 


2 


12 x 96 


256 


58 




13 


5 


3 


12 x 48 


512 


290 




14 


5 


4 


12 X 24 


1024 


754 




15 


5 


5 


12 X 12 


1024 


811 


15 


16 


6 


1 


6 x 192 


256 







17 


6 


2 


6 x 96 


512 


122 




IS 


6 


3 


6 x 48 


1024 


610 




19 


6 


4 


6 x 24 


2048 


1586 




20 


6 


5 


6 x 12 


4096 


3538 




21 


6 


6 


6x6 


4096 


3721 


18 


22 


7 


1 


3 x 192 


512 







23 


7 


2 


3 x 96 


1024 


252 




24 


7 


3 


3 x 48 


2048 


1260 




25 


7 


■1 


3 x 24 


4096 


3276 




26 


7 


5 


3 x 12 


8192 


7308 




27 


7 


6 


3x6 


16384 


15372 




28 


7 


7 


3x3 


16384 


15876 


21 



in region 20. This corresponds to scales of a few arcminutes 
(~ 6 — 12') which is expected since the characteristic signa- 
ture of the Kaiser-Stebbins effect shows at small scales. As 
an illustration, we also show in the middle panel the corre- 
sponding rotate d un correlated variables. These are obtained 
as explained in £3.2. 



To evaluate in a more quantitative way the discriminat- 
ing power of our non-Gaussianity test, we have obtained the 
Fisher discriminant function t as explained in § 3.3. This is 
shown in the bottom panel of Fig. ^. The solid and dashed 
lines correspond to the histogram of t for the Gaussian and 
non-Gaussian simulations respectively. As expected from the 
difference in the fourth cumulant spectrum, the curves are 
clearly separated. Indeed the power of the test to discrim- 
inate between both models at the 0.01 significance level is 
1. 

In order to assess the amount of Gaussian confusion that 
can be added to our test map and still detect the underlying 
non-Gaussianity, we increase the value of b until the mod- 
els cannot be separated. Fig. ^ shows the fourth cumulant 
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Table 2. Discriminating power of the non-Gaussianity test for different wavelets 



Tensor Mallat A trous 



Daubechies 4 Haar Daubcchics 4 Haar 







0° 


180° 




0° 


180° 






b = 





1.00 


1.00 


1.00 


1.00 


1.00 


1.00 


1.00 


b = 


1 


1.00 


1.00 


1.00 


1.00 


1.00 


1.00 


1.00 


b = 


2 


0.517 


0.798 


0.619 


0.974 


0.996 


0.995 


0.991 


b = 


3 


0.038 


0.091 


0.037 


0.165 


0.330 


0.252 


0.187 



b=1 




region number 



b=2 




b=3 




region number 



Figure 3. Fourth cumulant spectra and Fisher discriminant for b=l,2,3 using the tensor Daubechies 4 wavelet. Open circles/solid lines 
and solid squares/dashed lines correspond to the Gaussian and non-Gaussian models respectively. 



spectrum and Fisher discriminant function for b — 1,2,3. 
For non-Gaussian:Gaussian proportions of 1:1, the presence 
of the non-Gaussian signature is still clear at region 20, 
although the level of the detection has fallen appreciably. 
This translates into much closer histograms for the Fisher 
discriminant function but the power of the test at the 0.01 
significance level remains 1. When the level of Gaussian con- 



fusion is increased to b = 2, the fourth cumulant spectrum of 
both models shows no clear detections. Although some dif- 
ferences are seen in the estimated mean values of both popu- 
lations, the error bars clearly overlap. However, when all the 
information is taken into account simultaneously using the 
Fisher discriminant, the models can be partially separated. 
We find that the power of the test is p = 0.517. This means 
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x 2 

Figure 4. Histogram of \ 2 values as derived from the fourth 
cumulant spectra of the Gaussian (solid line) and non-Gaussian 
(dashed lines) simulations for 6 = 2 when using tensor Dau bcch ics 
4 wavelets. The x 2 values were obtained as explained in §p,2| 



that we are able to detect the non-Gaussian signal in ap- 
proximately half of our simulated test maps. For 6 = 3, the 
non-Gaussianity signal is completely obscured by the Gaus- 
sian component, and the test fails to discriminate between 
both models (p = 0.038). Table ^ summarises the values of 
p for the different cases considered. 

As pointed out in the previous section, the Fisher func- 
tion is the optimal linear discriminant in the sense of max- 
imising the power of the test. It is interesting to compare 
how the standard x 2 test performs in comparison with the 
previous results. To illustrate this point, we plot the his- 
togram for the x 2 values obtained from the Gaussian and 
non-Gaussian simulations for the case b = 2 (see Fig.g). We 
see that in this case the curves for both models completely 
overlap and the power of the test drops to p — 0.017 to be 
compared with p — 0.517 in the case of the Fisher discrimi- 
nant function. Indeed, we find that the Fisher discriminant 
function clearly outperforms the power of the \ 2 test in most 
cases. Only in those cases where the differences between the 
models are very large (for instance when b — 0), do both 
tests give comparable results. 

Another point that needs to be addressed is the orien- 
tation sensitivity of the test due to the assymetry of the 
wavelet basis functions. Because of the homogeneity and 
isotropy of the CMB temperature field, different orientations 
of the data should be statistically equivalent. However, if we 
look at one single realisation the results will vary when ro- 
tating the data. This is also what is happening in our case, 
since we are using the same underlying cosmic strings real- 
isation for all the non-Gaussian maps. Therefore, repeating 
the former analysis for a different orientation of our simu- 
lated maps will give us an insight on the importance of this 
effect when analysing a particular set of real data. Thus, 
we have performed the analysis for the same Gaussian and 
non-Gaussian maps rotated 180 degrees with respect to the 
original orientation. The discriminating power of the test for 
these cases are given in Table ^. These numbers are typical 
of the range of variation of the power of the test for differ- 
ent rotations. For illustration, we show the case for b — 1 in 
Fig. P We see that the single detection of non-Gaussianity at 
region 20 is larger than that of the original orientation (see 
middle panels of Fig.^J) . However, when all the information 
is taken into account, the overlapping between both mod- 



b=1,rotated 



it tdlH 




Figure 5. Fourth cumulant spectrum and histogram of Fisher 
discriminant for the Gaussian and non-Gaussian simulations (b = 
1) rotated by 180 degrees with respect to the original orientation. 
Tensor Daubechies 4 wavelets were used. 



els is quite similar. For larger levels of confusion (b = 2, 3), 
the Fisher discriminant function is sensitive even to small 
differences in the detection level of non-Gaussianity, leading 
to a better separation of the models for the rotated data. In 
particular, p = 0.798 for b = 2 versus p = 0.517 for the origi- 
nal orientation of the maps. Therefore, our non-Gaussianity 
test is sensitive to the orientation of the data. In order to 
deal with this problem, we need to use symmetric (or anti- 
symmetric) wavelet basis functions. For orthogonal wavelets, 
only the Haar wavelet fulfils this requirement. Thus, in the 
next subsection we repeat our analysis using this wavelet. 
A second possibility is to use a symmetric non-orthogonal 
wavelet, su ch a s the one constructed with the a trous algo- 
rithm (see §2.3). 



4.1.2 Haar 

As already pointed out, using the Haar wavelet solves the 
problem of the orientation sensitivity of our non-Gaussianity 
test. Indeed, with this wavelet basis, our test is invariant 
under rotations of 90, 180 and 270 degrees of the original 
data. Moreover, there are no edge effects when using the 
Haar wavelet, and therefore all the wavelet coefficients can 
be used in the analysis. Haar is also the simplest orthogonal 
wavelet to implement. In Fig. ^|we plot the fourth cumulant 
spectra and the histograms of the Fisher discriminant for 
b — 0,1,2,3. We see that the Gaussian and non-Gaussian 
models are successfully separated for 6 = and 6=1, with 
a discriminating power of 1 (see Table |^) . When the level of 
Gaussian confusion is increased the models start to overlap 
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Figure 6. Fourth cumulant spectrum and Fisher discriminant for different values of b using tensor Haar wavelets. 



(p — 0.619 for b — 2) and finally mix completely (p = 0.037 
for 6 = 3). These results are very similar to those obtained 
using Daubechies 4 (see Figs, g and We note that the 
discriminating power of the test for b — 2 in this case lies in 
between the values obtained for the two different rotations 



of the data when the analysis is performed with Daubechies 
4. 

Therefore the Haar wavelet gives similar results to the 
Daubechies 4 case, but it also has the property of rotational 



mvana.iice. 
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Table 3. Correspondence of region number with scale for Mallat 
algorithm. The scale is calculated as 1.5'(2 J— 3 ). The total number 
of coefficients and those used in the analysis of each region when 
using Daubechies 4 are given in the fifth and sixth columns. All 
the coefficients are used for the Haar wavelet. The corresponding 
region number for tensor wavelets of those regions common to 
both kind of wavelets are given in the last column. 



Table 4. Correspondence of region number with scale for the d 
trous algorithm. The scale is estimated as 1.5' (2 3 ) The last 
two columns give the total number and the number of coefficients 
used to calculate the fourth cumulant at each region 



region 


j 


detail 


~scale(') 


no.coeff. 


used (D4) 


Tensor 


1 


1 


V 


192 


4 







2 


1 


H 


192 


4 







3 


1 


D 


192 


4 





1 


4 


2 


V 


96 


16 







5 


2 


H 


96 


1(5 







6 


2 


D 


96 


1(5 





3 


7 


3 


V 


48 


61 


25 




8 


3 


H 


48 


64 


25 




9 


3 


D 


48 


64 


25 


(5 


10 


4 


V 


24 


256 


169 




11 


4 


H 


24 


256 


1(59 




12 


1 


D 


24 


256 


169 


10 


13 


5 


V 


12 


1024 


841 




11 


5 


H 


12 


1024 


841 




15 


5 


D 


12 


1024 


841 


15 


16 


6 


V 


6 


4096 


3721 




17 


6 


H 


6 


4096 


3721 




18 


6 


D 


6 


4096 


3721 


21 


19 


7 


V 


3 


16384 


15876 




20 


7 


H 


3 


16384 


15876 




21 


7 


D 


3 


16384 


15876 


28 



type of wavelet, there are three different kind of coefficients: 
vertical, horizontal and diagonal. Horizontal and vertical de- 
tails should be statistically equivalent, but we obtain differ- 
ent levels of detection because we are looking at a particular 
realisation (the same underlying non-Gaussian map), as it 
would happen for a set of real data. The numbering of the 
regions is as follows: regions 1,2,3 correspond to vertical, 
horizontal and diagonal details respectively for the largest 
scale (j — 1); regions 4,5,6 give the three details in the same 
order for the next scale (j = 2) and so on (see Table pi). 



4-2.1 Daubechies 4 

Figure ^ shows the fourth cumulant spectra for different 
non-Gaussian:Gaussian proportions (b = 0,1,2,3) as well 
as the histograms of the Fisher discriminant function for 
the Gaussian and non-Gaussian models. The values of the 
power of the test are summarised in Table ^. Again we see 
how the models begin to overlap as the level of the Gaussian 
component is increased. However, the test is more powerful 
at discriminating between the models than the Daubechies 
4 tensor case (see Figs. ^| and ^). There are a larger num- 
ber of clear detections which also have a higher amplitude. 
In particular, the highest detection is found at region 16, 
which corresponds to the vertical detail coefficients at scales 
of ~ 6'. The higher discriminating power of these wavelets 
becomes especially clear for the case 6 = 2, where we find a 



region 


3 


~scale(') 


no.coeff. 


used 


2 


2 


96 


65536 


16 


3 


3 


48 


65536 


17424 


■1 


4 


24 


65536 


38416 


5 


5 


12 


65536 


51984 


6 


(5 


6 


65536 


59536 


7 


7 


3 


65536 


63504 



4.2 Mallat algorithm 

In this section we repeat the previous analysis for 
Daubechies 4 and Haar but using 2-dimensional wavelets 
constructed as explained in §2.2. As already said, for this 



value of p = 0.974 versus p = 0.519 for the tensor wavelets. 
We note that some of the regions are common for both the 
tensor and Mallat wavelets (see Table |^). However, the high- 
est detections of non-Gaussianity do not appear in these 
regions. In particular, for the Mallat algorithm, the non- 
Gaussian detection in the horizontal and vertical details is 
usually larger than that of the diagonal wavelet coefficients 
for the same scale. This may be due to the type of non- 
Gaussianity present in our strings map, which contain sharp 
features that are more easily detected on the horizontal and 
vertical details. Analogously, for the tensor case, the largest 
detections are obtained in regions with slightly different val- 
ues of ji and j2, since they seem to identify better the fea- 
tures of the long strings. 

The corresponding values for a different orientation of 
the data (rotated 180 degrees with respect to the original 
maps) are given in Table ti. As for the tensor case, the single 
detections of non-Gaussianity can appreciably vary for dif- 
ferent rotations. However, when all the information is taken 
into account in the Fisher discriminant, they lead to similar 
results, at least for the cases of moderate levels of confusion. 



4.2.2 Haar 

In order to address the orientation sensitivity problem, we 
repeat the former analysis again using the Haar wavelet. The 
fourth cumulant spectra and Fisher discriminant function 
for 6 = 0, 1, 2, 3 are given in Fig. |^. Although the amplitude 
of the single non-Gaussian detections is different from the 
Daubechies 4 case (see Fig. [?]), the discriminating power 
of the test is actually very similar. In fact, the values of p 
for the Haar wavelet, lies between the values for the two 
considered orientations of the data. We find again that the 
Mallat wavelet works better discriminating between both 
models than the tensor construction when the same wavelet 
basis is used (see Table ^) . 



4.3 A trous algorithm 

In this section we compute the results for the same test but 



now based on the a trous algorithm. As explained in §2.3, 
this is a redundant non-orthogonal wavelet basis, that con- 
tains the same number of wavelet coefficients at each scale as 
the original image. Since we have chosen a symmetric mask, 
the a trous wavelet transform is invariant under rotations 
of the data of 90, 180 and 270 degrees with respect to the 
original orientation. Fig. ^ shows the results for the fourth 
cumulant spectra and the Fisher discriminant when the a 
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Figure 7. Fourth cumulant spectrum and Fisher discriminant for different values of b using Mallat Daubechies 4 wavelets. 



trous algorithm is used for b = 0, 1, 2, 3. The region number 
coincides with the j scale of the wavelet. Therefore, higher 
region numbers correspond to smaller scales (see Table |^). 
We see that, the highest single detection occurs at region 6, 
that corresponds to scales of ~ 6'. This coincides with the 
scale where the maximum detection was found for the former 



studied wavelets. As an illustration, Fig. [L0| also shows one 
of our non- Gaussian (top left panel) and Gaussian (top right 
panel) temperature maps for 6=1 and the corresponding 
wavelet coefficients (bottom panels) at scale j — 6, where 
the maximum detection is found. We see clear differences 
in the structure of the wavelet coefficients maps, which are 
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Figure 8. Fourth cumulant spectrum and Fisher discriminant for b = 0, 1,2, 3 for Mallat Haar wavelets. 



however not appreciated in the maps themselves. In fact, 
the wavelet transform succeeds on picking up the butterfly- 
like signature of the Kaiser-Stebbins effect which is seen in 
the wavelet coefficients of the non-Gaussian map. In spite of 
having only six different scales, that contain repeated infor- 
mation (since each region is comprised of the same number 



of wavelet coefficients as the original image) the separation 
of the models is better than that obtained with the tensor 
wavelets and comparable to the Mallat algorithm case (see 
Table §). 

In HJL, the case where Gaussian white noise is the main 
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Figure 9. The fourth cumulant spectra and the corresponding histograms of the Fisher discriminant are given for various values of b 
obtained when using the d trous algorithm. 



contaminant to the CMB map was also considered. This has 
been traditionally the greatest obstacle to the detection of 
non-Gaussianity in the CMB, and this is still the case for the 
4-yr COBE data. However, current and future experiments 
are expected to produce data with a much higher signal to 



noise ratio and therefore this case offers less interest. Due 
to this fact and for the sake of brevity, we do no include 
this case in the present work. Nonetheless, we have repeated 
our analysis for maps where the rms level of instrumental 
noise is equal to that of the convolved CMB map and find 
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that all the considered wavelets perform very well. In fact 
the two models are clearly discriminated being the power 
of the test p = 1. We also find that the highest detection 
move towards larger scales. This is expected since the noise 
is most contributing at small scales, therefore obscuring the 
non-Gaussian signal. 



5 DISCUSSION AND CONCLUSIONS 

We have investigated the performance of wavelet techniques 
to detect and characterise non-Gaussianity in the CMB, ex- 
tending the work of Hobson, Jones & Lasenby (1999). Our 
test map consists of a realisation of the Kaiser-Stebbins ef- 
fect due to cosmic strings of size 6.4 square degrees and pixel 
of 1.5 arcminutes. In addition, a Gaussian component with 
the same power spectrum as the original signal has been 
superposed with Gaussiamnon-Gaussian proportions 1 : b 
and the resulting signal convolved with a Gaussian beam 
of FWHM = 5 arcminutes. Gaussian white noise was also 
added with an rms equal to one-tenth that of the smoothed 
signal. A large number of non-Gaussian maps (10000) were 
generated, containing the same underlying Kaiser-Stebbins 
map and different realisations of the Gaussian component 



and noise. This allows us to account for the dispersion due 
to the Gaussian component and noise although the results 
will depend on the particular non-Gaussian underlying cos- 
mic string map that is being used. Ideally, different non- 
Gaussian components should have been used for each map 
but generating cosmic strings maps is very computationally 
expensive and a large number of them were not available to 
us. However, these maps are still a good test for the method 
presented and, even if the particular level of discrimination 
may vary from realisation to realisation, we expect the gen- 
eral conclusions of this work to remain valid. The same num- 
ber of Gaussian simulations were generated, simply by ran- 
domising the phases of the non-Gaussian maps. This method 
ensures that Gaussian and non- Gaussian realisations share 
the same power spectrum and thus any difference found be- 
tween both models will be due to the non-Gaussian character 
of the signal and not to differences in the power spectrum. 

The non-Gaussian test is as follows. First, we have ob- 
tained the wavelet transform of our non-Gaussian and Gaus- 
sian maps. We have then computed the fourth cumulant at 
each wavelet scale for each of the maps. All this information 
has been combined by computing the Fisher discriminant 
function for each simulated map and, finally, the distribu- 
tion of this quantity has been compared for the Gaussian 
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and non- Gaussian models. Indeed, the choice of a good test 
statistic that combines all the information from the different 
wavelet scales is crucial. The Fisher discriminant function is 
the optimal linear function of the measured variables in the 
sense of maximum power. In fact, we have found that, for 
most of our cases, the Fisher discriminant function clearly 
outperforms a generalised x 2 test. 

We have also compared the performance of different 
wavelet transforms. Wavelets constructed using the Mallat 
algorithm have been shown to perform better than the 2D- 
tensor wavelets (used in the original work of HJL). In partic- 
ular, using the Mallat algorithm, the stringy and inflation- 
ary models can be discriminated at a very high confidence 
level even when the Gaussian confusion has a rms twice as 
large as the one of the cosmic strings map (6 = 2, p — 0.995 
for Haar). However, the discriminating power of the test is 
clearly lower for the same case when tensor wavelets are used 
(p = 0.619). In addition, we have tried different orthogonal 
bases for both decomposition schemes and find very simi- 
lar results, providing all the information is combined into 
the Fisher discriminant function. In general, we have found 
that simpler wavelet bases perform slightly better than more 
complex ones. Thus we have used the Daubechies 4 and Haar 
wavelets for our work. 

We have also investigated the performance of the a trous 
algorithm, a non-orthogonal wavelet transform. In spite of 
being redundant (each scale contains as many coefficients as 
the original image), its discriminating power is very similar 
to that of the Mallat algorithm. In particular, we find p = 
0.991 for 6 = 2. 

We have also addressed the problem of orthogonal 
wavelet transforms being sensitive to the orientation of 
the data, which is a consequence of the assymetry of the 
wavelet basis functions. To avoid this problem, we have used 
the Haar wavelet functions, which are antisymmetric. This 
means that the non-Gaussianity test is invariant under ro- 
tations of the data when the Haar wavelet is used. Another 
possibility is using a symmetric non-orthogonal wavelet such 
as the a trous algorithm. We have also presented results for 
Daubechies 4 for two different orientations of the data to 
give an insight of the sensitivity of the test when the origi- 
nal map is rotated. 

In addition, wavelet techniques are useful to charac- 
terise the angular scale at which the non-Gaussian signal 
occurs. We have found the maximum non-Gaussian signal 
at scales ~ 6 arcminutes for the different considered trans- 
forms. Finally, we note that our test is not making use of the 
localisation property of wavelets. However, we can clearly 
see in the wavelet coefficients map of Fig. llOl individual pat- 
terns of the Kaiser- Stebbins effect. This shows that further 
wavelet algorithms can be investigated in order not only to 
detect the non-Gaussian signature of cosmic strings globally 
but also to extract individual non-Gaussian structures from 
the map. 
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